%% prepare plot3d file for gridgen mesh

x1=-1000;x2=1000;
y1=-1000;y2=1000;
z2=500; zm=1000;

nx=20;ny=20;
[x,y]=meshgrid(linspace(x1,x2,nx),linspace(y1,y2,nx));

z=exp(-(sqrt(x.^2+y.^2))/abs(x2));
z=-z*zm;


fid=fopen('surf.net','w');
fprintf(fid,'%d\n',6);
fprintf(fid,'%10d %10d %10d ',[size(z),1]); %surface
fprintf(fid,'%10d %10d %10d ',[nx,2,1]);    
fprintf(fid,'%10d %10d %10d ',[nx,2,1]);
fprintf(fid,'%10d %10d %10d ',[ny,2,1]);
fprintf(fid,'%10d %10d %10d ',[ny,2,1]);
fprintf(fid,'%10d %10d %10d\n',[nx,ny,1]);
%up
fprintf(fid,'%15f %15f %15f %15f \n',x);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',y);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',z);
%front
index=1:ny:nx*ny;
fprintf(fid,'%15f %15f %15f %15f \n',[x(index),x(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[y(index),y(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[z(index),z(index)*0+z2]);
%back
index=ny:ny:nx*ny;
fprintf(fid,'%15f %15f %15f %15f \n',[x(index),x(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[y(index),y(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[z(index),z(index)*0+z2]);
%left
index=1:ny;
fprintf(fid,'%15f %15f %15f %15f \n',[x(index),x(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[y(index),y(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[z(index),z(index)*0+z2]);
%right
index=(nx-1)*ny+(1:ny);
fprintf(fid,'%15f %15f %15f %15f \n',[x(index),x(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[y(index),y(index)]);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',[z(index),z(index)*0+z2]);
%bottom
fprintf(fid,'%15f %15f %15f %15f \n',x);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',y);
fprintf(fid,'\n',1);
fprintf(fid,'%15f %15f %15f %15f \n',z*0+z2);

fclose(fid);












